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We develop convergent variational perturbation theory for quantum statistical density matrices. 
The theory is applicable to polynomial as well as nonpolynomial interactions. Illustrating the power 
of the theory, we calculate the temperature-dependent density of a particle in a double-well and of 

bo", 

^\ ' the electron in a hydrogen atom. 
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(N ! I. INTRODUCTION 



' Variational perturbation theory [jlj|2|] transforms divergent perturbation expansions into convergent ones. The 
convergence extends to infinitely strong couplings 0], a property which has recently been used to derive critical 



exponents in field theory without renormalization group methods 0J5|. The theory has first been developed in 

00 ■ . . 

• quantum mechanics for the path integral representation of the free energy of the anharmonic oscillator |3| and the 

>: 

i 



hydrogen atom Local quantities such as quantum statistical density matrices have been treated so far only to 

lowest-order for the anharmonic oscillator and the hydrogen atom H||. There has also been a related first-order 



& . 

^ , treatment in classical phase space |Kj] for systems with dissipation . 

a -1 

The purpose of this paper is to develop a systematic convergent variational perturbation theory for the path integral 

• i-H , 

representation of density matrices of a point particle moving in a polynomial or nonpolynomial potential. As a first 
application we calculate the particle density in a double-well and then the electron density in a hydrogen atom. 



II. GENERAL FEATURES 



Variational perturbation theory approximates a quantum statistical system by perturbation expansions around 
harmonic oscillators with trial frequencies which are optimized differently for each order of the expansions. When 
dealing with the free energy, it is essential to give a special treatment to the fluctuations of the path average x = 
(ksT/Ti) j^ ksT dr x(t), since this performs violent fluctuations at high temperatures T. These cannot be treated 
by any expansion, unless the potential is close to harmonic. The effect of these fluctuations may, however, easily be 



1 



calculated at the end by a single numerical fluctuation integral. For this reason, variational perturbation expansions 
are performed for each position xo of the path average separately, yielding an Nth order approximation Wn(xq) to 
the local free energy T4ff, c i(xo), called the effective classical potential Il2]. The name indicates that one may obtain 
the full quantum partition function Z from this object by a simple integral over xo just as in classical statistics, 

/+oo 1 
X ° exp{-V/ eff , cl (x )/fcBT}. (2.1) 

-°° ^2Trh 2 /Mk B T 

Having calculated Wn(xo), we obtain the TVth-order approximation to the partition function 

+00 

Z N = j e -W N (,o)/k B T (22) 

i ^2irh 2 /Mk B T 

The separate treatment of the path average is important to ensure a fast convergence at larger temperatures. In the 
high-temperature limit, Wn(xq) converges against the initial potential V(xq) for any order N. 

Before embarking upon the theory, it is useful to visualize some characteristic properties of path fluctuations. 
Consider the euclidean path integral over all periodic paths x(t), with x(0) — xfi/ksT), for a harmonic oscillator 
with minimum at x TO , where the action is 

An, x Jx] = J oIt |-Mx 2 (r) + -MH 2 [x{t) - x m } 2 j . (2.3) 

Its partition function is 

z ^ = j Vx «p{-^W/»> - 2 ^m,2k B T (2 - 4) 

and the correlation functions of local quantities 0\ (x) , O2 (x) , . . . are given by the expectation values 

(0 1 (x(r))0 2 (x(r))--- ) n > x ™ = -i— /dx0 1 (x(t 1 ))0 2 (x(t 2 )) • ■ -exp{-.An,*Jx]/7i} . (2.5) 



The particle distribution of the oscillator is given by 



(x) = (S(x - x(t)) ) n <*™ = -J== exp 



x m ) 



(2.6) 



which is a Gaussian distribution of width 



a„ = coth , (2.7) 

H 2MO 2fc B T' V ' 

the subscript indicating that we are dealing with a harmonic oscillator. At zero temperature, this is equal to the 
square of the ground-state wave function of the harmonic oscillator, whose width is 
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FIG. 1. Temperature dependence of fluctuation widths of any point x(r) on the path in a harmonic oscillator (I 2 is a square 
length in units of %/MQ). The quantity (dashed) is the quantum mechanical width, whereas a 2 , Q (dash-dotted) shares the 
width after separating out the fluctuations around the path average xq. The quantity Og cl (long-dashed) is the width of the 
classical distribution, and fe^ (solid curve) is the fluctuation width at fixed ends which is relevant for the calculation of the 
density matrix by variational perturbation theory. 



In the limit fi — > 0, we obtain from (2.5), ( |2.6| ) the classical distribution 

{X ) 



Pnci(x) 



1 



: exp 



24 cl 



(2.9) 



with 



l Hcl 



Mfl 2 



(2.10) 



The linear growth of this classical width is the origin of the famous Dulong-Petit law for the specific heat of a harmonic 
system. The classical fluctuations are governed by the integral over the Boltzmann factor 



-Mn 2 (x-x m ) 2 /2k B T 



(2.11) 



in the classical partition function 



3 



ZhcI — 



dx 



2Trh 2 /Mk B T 



-Mn 2 (x-x m ) 2 /2k B T 



(2.12) 



From this we obtain the classical distribution (2.9) as the expectation value 

+00 



PhciO) = (S(x- x))ci Xm = Z cl 



dx 



2irh 2 /Mk B T 



: 8{x — x)e 



-\ -MQ, 2 (x-x m ) 2 /2k B T _ 



y/Tl 



: exp 



2ai 



(2.13) 



Variational perturbation theory avoids the divergence of the harmonic width at high temperatures ( 2.10 ) by the 
separate treatment of the fluctuations of the path average x, as explained above. The average is fixed at some value 
Xq with the help of a delta function S(x — xq). For each xq we introduce local expectation values 

( S(x - x Q ) Oi(x(n)) 2 (x{t 2 )) ■ ■ ■ ) n '*™ 



Sl,x„ 



(2.14) 



(5(x - X ) ) n ' x ™ 

The original quantum statistical distribution of the harmonic oscillator ( |2.6| ) collects fluctuations of x — xq and those 
around xq, and can therefore written as a convolution 



+ 00 



dXQ Px {x - X ) Pcl{x ), 



(2.15) 



over the classical distribution (2.9) and the local one 



P xo (x) = (6(x - x(r)))l 



exp 



(X - Xq) 2 
2< 



(2.16) 



Such a convolution of Gaussian distributions (2.15) leads to another Gaussian distribution with added widths, so that 



the width of the local distribution is given by the difference 



— an — a^, = — cotii — 

Xo H cl 2MQ \ 2k R T hVL 



nn 2k R T 



(2.17) 



which starts out at a finite value for T = 0, and goes to zero for T — > 00, 

o HQ 
lim ~- - 



(2.18) 



x " \2k B T' 

The latter property suppresses all fluctuations around x and guarantees that limT->oo Wn{xq) = V(xq) for all N (see 
Fig. 0). 

With this separation of the path average, the partition function 



Z = <j>Vx exp{-A[x]/H} 



(2.19) 



for the general particle action 



A[x] = I 
Jo 



h/k B T 



-Mx 2 (t) + V(x(t)) 



(2.20) 



possesses the effective classical representation (pjl) with the effective classical potential 



Veff.dfxo) = -k B T In 



2tt?i 



Vx5{xq — x) exp {— „4[x]/7i} 



(2.21) 



Mk B T , 

In variational perturbation theory, this is expanded perturbatively around an xo-dependent harmonic system with 
trial frequency Q(xq), whose optimization leads to the approximation Wn(xq) for V e s^ c i(xo)- 



III. DENSITY MATRIX OF HARMONIC OSCILLATOR 



How can this method be extended to density matrices? Their path integral representation is 



p(x b ,x a ) = —p(x b ,x a ) 



(3.1) 



where p(x b ,x a ) is the path integral 



p(x b ,x a ) = J Vx exp {— A[x}/h} (3-2) 

(x a ,0y* (x b ,h/k B T) 

over all paths with the fixed endpoints x(0) = x a and x{h/k B T) = x b . The partition function is found from the trace 
of p(x b ,x a ): 

+ 00 

Z = J dx x). (3-3) 

— oo 

For a harmonic oscillator centered at x m ( |2.3| ), the path integral (3.2) can be easily done with the result [^) 
Po (x b ,x a ) = 



Am 



■ exp 



where 



2nhsmhhn/k B T { 2hsmhhfl/k B T 



x(t) = x(t) — X Tl 



[(xl + x 2 .) coshhtl/ k B T - 2x b x a ] 



(3.4) 



(3.5) 



At fixed endpoints x b ,x a , the quantum mechanical correlation functions are 

Vx0 1 (x(t 1 ))0 2 (x(t 2 )) ■ ■ • exp {-Aw MAI 



(OiWnDO,^))...)^;^ 1 



Po {x b ,X a ) 



{x a ,0)~->{x b ,h/k B T) 



(3.6) 



and the distribution function is given by 

Ph (x,t) = (8(x-x(t)))°1*£ = 1 exp 

v / 2tt6 h (t) 

The classical path of a particle in a harmonic potential is 



(x - x c i) 2 
26^(r) 



(3.7) 



i h sinhfir + i a sinhil(?i/fc B r - t) 
* cl(T) = sinhMV^T (3 - 8) 



and the time-dependent width &h( t ) * s f° un d to be 



, _ h ( TiCl cosh[fi(2r - h/k B T)] 



6h(t) = im ( COth sinhMVfc B r * • (3 - 9) 



Since the euclidean time r lies in the interval < r < h/ksT, the width (3.9) is bounded by 



thus remaining finite at all temperatures. The temporal average of ( |3.9| ) is 

, 2 /c b t /- R / fe B T n f nn k B T\ 

b ^nrj dTh - {T) = ^{ coth k^-in)- (3 - n) 



Just as a^ , this goes to zero for T — > oo with an asymptotic behaviour hfl/6ksT, which is twice as big as that of a^ o 
(see Fig. |). 

IV. VARIATIONAL PERTURBATION THEORY FOR DENSITY MATRICES 



To obtain a variational approximation for the density matrix, it is useful to separate the general action ( 2.20| ) into 



a trial one for which the euclidean propagator is known, and a remainder containing the original potential. If we were 
to proceed in complete analogy with the treatment of the partition function, we would expand the euclidean path 
integral around a trial harmonic one with fixed end points xt,, x a and a fixed path average Xo, and with a trial frequency 
Q(xb, x a ;xo). The result would be an effective classical potential Wn{xb, x a ]Xo) to be optimized in £l(xb, x a ]Xo). After 
that we would have to perform a final integral in xq over the Boltzmann factor exp[— Wjyfab, x a ; Xaj/ksT]. 



G 




k B T/hn 



FIG. 2. Temperature-dependence of the width of fluctuations around the path average xq = x at fixed ends. For comparison 
we also show the width a^ of Fig. [j]. The vertical axis gives these square lengths in units of h/Mil again. 

But, because of the fmiteness of the fluctuation width 6^ at all temperatures which is similar to that of a 2 , , the 
special treatment of x = xq becomes superfluous for paths with fixed endpoints Xb,x a - While the separation of xq was 
necessary to deal with the diverging fluctuation width of the path average x, paths with fixed ends have fluctuations 
of the path average which are governed by the distribution 

1 f 1 



Px (x b ,x a ;x ) = (5(x-x)) 



: exp ■ 



2&? 

Xq 



\, s 2k B T , nn 



(4.1) 



with the width 



2 = k B T 
Xo MO 2 



2k B T , Ml 
1 — — tanh ■ 



(4.2) 



Ml 2k B T 

which goes to zero for both limits T — > and T — > oo (see Fig. ^). At each euclidean time, x(t) fluctuates narrowly 
around the classical path x c \(t) connecting Xb and x a . This is the reason why we may treat the fluctuations of x — x$ 
by variational perturbation theory, just as the other fluctuations. As a remnant of the extra treatment of xq we must, 
however, perform the initial perturbation expansion around the minimum of the effective classical potential which will 
lie at some point x m determined by the endpoints Xb,x a , and by the minimum of the potential V(x). Thus we shall 
use the euclidean path integral for the density matrix of the harmonic oscillator centered at x m as the trial system 
around which to perform the variational perturbation theory, treating the fluctuations of Xq around x m on the same 
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footing as the remaining fluctuations. The position x m of the minimum is a function x m = x m (x b , x a ), and has to be 
optimized with respect to the trial frequency, which itself is a function f2 = Q(x b ,x a ) to be optimized. 



Hence we start by decomposing the action (2.20) as 



A[x] = An, Xm [x] + Aint[x] 



(4.3) 



with an interaction 



A int [x(T)} 



ti/) 



dr V5 n t (a; (t)), 



(4.4) 



where the interaction potential is the difference between the original one V(x) and the inserted displaced harmonic 
oscillator: 



V int (x(r)) = V(x(r)) - ^Mn 2 [x(r) - x m } 2 . 



(4.5) 



For brevity, we have introduced the inverse temperature in natural units (3 = 1/ksT in (4.4). Now we evaluate 



the path integral for the euclidean propagator (3.2) by treating the interaction (4.4) as a perturbation, leading to a 
moment expansion 



p(x b ,x a ) = p ,Xm (x b ,x a ) 



(4.6) 



with expectation values defined in (3.6). The zeroth order consists of the harmonic contribution (3.4) and higher 



orders contain harmonic averages of the interaction (4.4). The correlation functions in (4.6) can be decomposed into 
connected ones by going over to cumulants, yielding 



p(x b ,x a ) = p ' Xm (x b ,x a )exp 



(4.7) 



where the first cumulants are defined as usual 

(OiWnWtwHOiWn)))^, 
<Oi(x(r 1 ))0 2 (x(r 2 )))^;: c =(0 1 (x(t 1 ))0 2 (x(t 2 )))% x Z - <Oi(z(n (0 2 (x(r 2 )))^r , 



(4.8) 



The series ( |4. 7| ) is truncated after the iV-th term, resulting in the N-th order approximant for the quantum statistical 
density matrix 



p N (x b ,x a )=p ' (x b ,x a )exp 



N 



V I A? \x\) Q ' Xm 

n \ft n \^ l intL' t J lx b ,x a , 

71=1 



(4.9) 
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which explicitly depends on both variational parameters Q and x m . 

In analogy to classical statistics, where the Boltzmann distribution in configuration space is controlled by the 
classical potential V(x) according to 



1 M 

5— exp \—(3V(x)] , 

2tt?i 2 I3 



(4.10) 



we now introduce a new type of effective classical potential V e s tC \(x a ,x b ) which governs the unnormalized density 
matrix 



p{x b ,x a ) 



I M 

—2- exp [-PV cS . c i(x b , Xa)] 
2wh p 



(4.11) 



Its Nth. order approximation is obtained from (3.4), (4.9), and ( 4.11 ) voa the cumulant expansion 



W N ' Xm (x b ,x a ) = -^ln 



1 smhhfiVt 



Mfl 



(-1)" 



2(3 hpfl 2h/3smhh/3Q 



{{x 2 b + x 2 a ) coshh(3Q - 2x h Xa) - - y~] 

a *■ — ' nln 

n—l 



t Ant [ X ] )x b ,x a 



(4.12) 



which is optimized for each set of endpoints x b ,x a in the variational parameters fl 2 and x m , the result being denoted 
by WN(x b ,x a ). The optimal values fl 2 (x a ,x b ) and x m (x a , x b ) are determined from the extremality conditions 



dW£' Xm (X b , Xa)_ ± Q dWj?' X ">{ Xb ,X a ) ^ Q 



dn 2 



dx r , 



(4.13) 



The solutions are denoted by £l 2N , xl^, both being functions of x b , x a . If no extrema are found, one has to look for the 



flattest region of the function ( 4.12 ), where the lowest higher-order derivative disappears. Eventually the normalized 
density matrix is obtained from 



n 2JV x N 

p N (x b ,X a ) = p N ' m (x b ,X a ), 



(4.14) 



where 



Zn — 



<r> 2 " T-™ 

j — ' m ( \ 



(4.15) 



In principle, one could also optimize the entire ratio (4.14), but this would be harder to do in practice. Moreover, the 
optimization of the unnormalized density matrix is the only option, if the normalization diverges due to singularities 



of the potential. This will be seen in Sect. VIII B 
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V. SMEARING FORMULA FOR DENSITY MATRICES 



In order to calculate the connected correlation functions in the variational perturbation expansion (4.9), we must 



find efficient formulas for evaluating expectation values (pjj) of any power of the interaction (4.4) 



~ Q,x m i \ 
Pq \ X bi x a) ; r 

X a ,0 



n/3 



dn Vi nt (x(ri) + x m ) 



exp \ --A<n,x m [x + x m ] 



(5.1) 



This can be done by an extension of the smearing formalism, developed in Ref. jTj. We rewrite the interaction 
potentials as 



+ OC 



+ OC 



\ i„i (•'•(>/ )+•'-,„) = / '/:/ r ilj( (:./+•(-,„ ) / ^ exp{iA;z;} exp 



ft/3 



dTi\i8(r — ti)x(t) 



(5.2) 



and introduce a current 



J(r) = J2 iKhS{T 



(5.3) 



1=1 



so that (5.1) becomes 



(Ant [ x ] )x' h ,x a 



Po 



^— -TT / dri I dziV in t(zi + x min ) / — exp{iA ; z;} 

™{Xb, X a ) J=1 ^0 J-oo J-oo ^T 1 " 



K Q > X ™[J]. (5.4) 



The kernel K n ' Xm [J] represents the generating functional for all correlation functions of the displaced harmonic 
oscillator 



.n,x m r 7i _ / us — J / dr 



iT 2 ' x '"[J] = / Pi exp 



o 



yj 2 (r) + ^Mfi 2 i 2 (T) + J(t) x(t) 



(5.5) 



For zero current J, this generating functional reduces the euclidean harmonic propagator (3.4) 



K il ^[J = 0}=p o n ' x -(x b> x a ). 



(5.6) 



For nonzero J, the solution of the functional integral ( p.5|) is given by 

-ft/3 -, /-ft/3 /-ft/3 



^ S2 ' Xm W =Po°" m (^^a)exp 



■~J o drJ{r)x ci {r) + -^ ^ drjf dr' J(r) G n (r, r') J(r') 



where x c i(r) denotes the classical path (3.8) and G (t, t') the harmonic Green function 

H cosh^(|r - r'| - h/3) - cosh ft (t + r' - h/3) 



G n (r,r') 



2Mft 



sinh ?i/3r2 



(5.7) 



(5.8) 



The expression (5.7) can be simplified by using the explicit expression ( p.3| ) for the current J. This leads to a 
generating functional 
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K a ^[J] = p n ' x ™(x b ,x a ) exp (^A T x cl - i A t Ga| , 



(5.9) 



where we have introduced the n-dimensional vectors A = (Ai, . . . , A n ) T , x c i = (x c i(ti), . . . , x c i(t„)) t with the super- 
script T denoting transposition, and the symmetric n x n-matrix G whose elements are Gki = G n (Tk,Ti). Inserting 



( p.9|) into (5.4), and performing the integrals with respect to Ai, . . . , A n , we obtain the n-th order smearing formula 
for the density matrix 

r hf3 



( -^int [ x ] 



n 



+00 



dri / dzi V int {zi + x m ) 



1=1 L" u " _0 ° 

1 / 1 ^ 

x — = exp < — - > 

V(2T) n dotG P \ 2^ 



2fc - ^ci(Tfe)] Gj., 1 [zi - x c \(n)] 



(5.10) 



The integrand contains an n-dimensional Gaussian distribution describing both thermal and quantum fluctuations 



around the harmonic classical path x c \{t) of Eq. (3.8) in a trial oscillator centered at x m , whose width is governed 



by the Green functions (5 



For closed paths with coinciding endpoints (xb — x a ), formula (5.10) leads to the n-th order smearing formula for 
particle densities 



p(xa) = ^p(x a ,x a ) = — j> VxS(x(t = 0) - x a ) OXp{~ A[x] / H} , 



(5.11) 



which can be written as 



( A^tN )x'*Z = nx , — 7 TI / dn 1 dzi v ^( z i + x ™) 

Pn' (Xa) , , JO J -00 



v / (2n) n+1 det a 2 



ex P ( - 9 HZ z k a u z i ( 5 - 12 ) 



k,l=0 



\i.x m i \ 1 1 
Po \ x a) i =1 J» 

with zq = x a - Here a denotes a symmetric (n + 1) x (n + l)-matrix whose elements a|j = a 2 (rfc, 77) are obtained from 
the harmonic Green function for periodic paths G n ' p (r, r') as (see Chapters 3 and 5 in M) 



a 2( T T >) = A G ^Pf T t >) = h co S hn(\T- T '\-hp/2) 
1 ' ' M y ' ' 2MQ. sinhTi^fi/2 



(5.13) 



The diagonal elements a 2 = a(r, t) represent the fluctuation width (2.7) which behaves in the classical limit like (2.1C) 
and at zero temperature like ( [2.8] ). 

Both smearing formulas ( |5.1C| ) and ( 5.12j ) allow in principle to determine all harmonic expectation values for the 
variational perturbation theory of density matrices in terms of ordinary Gaussian integrals. Unfortunately, in many 
applications containing nonpolynomial potentials, it is impossible to solve neither the spatial nor the temporal integrals 
analytically. This circumstance drastically increases the numerical effort in higher-order calculations. 
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VI. FIRST-ORDER VARIATIONAL RESULTS 



The first-order variational approximation gives usually a reasonable estimate for any desired quantity. Let us 
investigate the classical and the quantum mechanical limit of this approximation. To facilitate the discussion, we 



first derive a new representation for the first-order smearing formula (5.12) which allows a direct evaluation of the 
imaginary time integral. The resulting expression will depend only on temperature, whose low- and high-temperature 
limits can easily be extracted. 



A. Alternative Formula for First-Order Smearing 

For simplicity, we restrict ourselves to the case of particle densities and allow only symmetric potentials V(x) 
centered at the origin. If V(x) has only one minimum at the origin, then also x m will be zero. If V(x) has several 
symmetric minima, then x m goes to zero only at sufficiently high temperatures (see Ref. Q). 

To first order, the smearing formula ( |5.12 ) reads 



{Ai » t[x]) *^ = tfW)J dT J i Vintiz) vogo - < exp ^ ~* ~ ' ^ r - r "™ } ' (fU) 

— oo 

so that Mehler's summation formula 



hp +oo 

1 f f dz 1 f 1 (z 2 + x 2 )aoo — 2zx a dQi 



i r {x 2 + x i2 ){\ + b 2 )-Axx'b\ r i. 2 , 2 ,l^ 6" 

exp < -- 9 ,7 f = exp <^ — (x 2 + x' 2 ) } > H n (x)H n (x') 

V^T^Ij 2 \ 2(1 -b 2 ) J V \ 2 V ') 4-! 2"n! V ; V ' 



(6.2) 



n=0 

leads to an expansion in terms of Hcrmite polynomials H n (x) , whose temperature dependence stems from the diagonal 



elements of the harmonic Green function (5.13): 



Here the dimensionless functions Co are defined by 



h/J 



,2 







Inserting (5.13) and performing the integral over r, we obtain 

n 

E/n\ smhhpil{n/2 — k) 
\k) MOfn/2 - k) ' ( ' 5) 



(„) i \ ^ /n\ sinh h(3Q(n/2 — k) 

13 ~ 2" cosh™ /2f^ Q \k) h(3n(n /2-k) 



At high temperatures, these functions of j3 go all to unity, 
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whereas at zero temperature: 



lim C 



(n) 




lim Cf ] = 1, 



71 = 0, 
71 > 0. 



(6.6) 



(6.7) 




6.0 7i/3fi 8.0 



FIG. 3. Temperature-dependence of the first 9 functions Ci"', where j3 = I/HbT 



According to ( 4.1 2| ) , the first-order approximation to the new effective potential ( 4.1 2] ) is given by 



I , sinh^/30 Mfl , , hl3VL , N 
'■• +— ^tanh^ + ^K) 



M = -5 In 



2/3 ft/3fi 



with the smeared interaction potential 



(6.8) 



(6.9) 



It is instructive to discuss separately the limits f3 — ► and /3 — * 00 of dominating thermal and quantum fluctuations, 
respectively. 

B. Classical Limit of Effective Classical Potential 



In the classical limit /3 — ► 0, the first-order effective classical potential (p.8|) reduces to 



13 



Wf' C '(Za) = 7;Mn 2 X 2 a + lim Vg(x a ). 

2 p->o 



(6.10) 



The second term is determined by inserting the high-temperature limit of the fluctuation width ( 2.10 ) and of the 



polynomials (|6.6| ) into the expansion (3.3), leading to 

lim V$(x a ) = limX: ^H n (V^i.) / - ***(*) e-^^g, (^MO^*) • (6-11) 



fib 

Then we make use of the completeness relation for Hermite polynomials 



1 °° 1 

V 71=0 



(6.12) 



which may be derived from Mehler's summation formula (6.2) in the limit 6 — ► 1 , to reduce the smeared interaction 
potential V^(x a ) to the pure interaction potential (4.5): 



\imV$(x a ) = V int (x a ). 



(6.13) 



Recalling (4.5) we see that the first-order effective classical potential ( 6.10 ) approaches the classical one: 



lim W?' c \x a ) = V(x a ). 

13^0 



(6.14) 



This is a consequence of the vanishing fluctuation width b 2 ^ of the paths around the classical orbits. This property is 



universal to all higher-order approximations to the effective classical potential (4.12). Thus all corrections terms with 
ri > 1 must disappear in the limit [3 — > 0, 



H n=2 



(6.15) 



C. Zero- Temperature Limit 



At low temperatures, the first-order effective classical potential (|6.8|) becomes 



(6.16) 



The zero-temperature limit of the smeared potential in the second term defined in (3.9) follows from Eq. (|6.3D by 



taking into account the limiting procedure for the polynomials Cg in (|3.7|) and for the fluctuation width a 2 (pjl). 



Thus we obtain with Hq(x) = 1 and the inverse length k = y/MQ/h: 



+ OC , 

r K 2 

lim V^fia) = / dz\ —H (Kz) 2 exp{-K 2 z 2 }V int (z). 

I^oo J V 7T 



(6.17) 
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Introducing the harmonic eigenvalues 



e? = nnin 



(6.18) 



and the harmonic eigenfunctions 



, f 2\ V4 



we can reexpress the zero-temperature limit of the first-order effective classical potential flj.16 ) with (6.17) by 



(6.19) 



W^ m (x a ) = E$ + (^\V int \^). 



(6.20) 



This is recognized as the first-order harmonic Rayleigh-Schrodinger perturbative result for the ground state energy. 
For the discussion of the quantum mechanical limit of the first-order normalized density, 



f*™ dx a (x a ) exp { - i ( A int [x] )" a , Xa } 



(6.21) 



we proceed as follows. First we expand ( 3.21 ) up to first order in the interaction, leading to 

+00 



P\{x a ) = Po(x a ) 



l-jj (Ant[x] 



(6.22) 



Inserting (|3.5|) and (6.3) into the third term in (6.22), and assuming VL not to depend explicitly on x a , the x a -integral 



reduces to the orthonormality relation for Hcrmitc polynomials 

+00 

dx a H n (x a )H (x a )e~ x " = 5 n0 , 



1 



2 n nk/7r 



(6.23) 



so that the third term in (6.22) eventually becomes 



+00 



+00 



dx a po(x a ) (Aiat[x]) 



(3 I dz \ — Vi nt (z) exp{— k z } H (kz) 



(6.24) 



But this is just the n = -term of (6.3) with an opposite sign, thus cancelling the zeroth component of the second 
term in ( |6.22| ) , which would have been divergent for — ► 00. 

The resulting expression for the first-order normalized density is 



Pl{Xa) = Po(Xa) 



„ +OO . 

1 - l&l Cfl" 5 H n(^a) / dz J—V int (z) exp(- K V) H u {kz) 



(6.25) 



The zero-temperature limit of is from (6.7) and (6.18) 



lim PC^ = 1 



(6.26) 
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so that we obtain from (6.25) the limit 



Pi{x a ) = Poixa) 



11 f / K 

1_2 L oMcfl ^n H n( KX a) / dz d — V int (z) exp{- k 2 z 2 }H n (Kz) H (kz) 

n=l U ' « J V 7r 



(6.27) 



Taking into account the harmonic eigenfunctions (6.19), we can rewrite (6.27) as 



= l^0(^)| 2 = [^(^a)] 2 - 2^(».) E ^(X^^^MI 



n>0 



(6.28) 



which is just equivalent to the harmonic first-order Rayleigh-Schrodinger result for particle densities. 

Summarizing the results of this section, we have shown that our method has properly reproduced the high- and 



low-temperature limits. Because of relation (6.28), the variational approach for particle densities can be used to 
determine approximately the ground state wave function tpo(x a ) for the system of interest. 



VII. SMEARING FORMULA IN HIGHER SPATIAL DIMENSIONS 

Most physical systems possess many degrees of freedom. This requires an extension of our method to higher 
spatial dimensions. In general, we must consider anisotropic harmonic trial systems, in which the previous variational 
parameter J7 2 becomes aflx D-matrix with (i, v = 1,2, .... D. 

A. Isotropic Approximation 



An isotropic trial ansatz 



n 2 - n 2 s 



(7.1) 



can give rough initial estimates for the properties of the system. In this case, the n-th order smearing formula ( 5.12 ) 
generalizes directly to 



int I J / r a ,r a 



1 ™ 



tif) 



dn I d D zi V in t (zz) 



V(27r)™ +1 det( 



exp 



1 " 

9 Zk a k. 



kl^l 



k.l=0 



(7.2) 



with the D-dimensional vectors z; = {zu, z^i, ■ ■ . , 2dz) t '■ Note, that greek labels fj,, v, . . . = 1, 2, . . . , D specify spatial 
indices and latin labels k,l, . . . = 0, 1, 2, . . . , n refer to the different imaginary times. The vector z denotes r a , the 
matrix a 2 is the same as in Section [V|. The harmonic density reads 



27ra oo 



-D 



exp 



1 ° 
2^2~ X! 



(7.3) 
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B. Anisotropic Approximation 



In the discussion of the anisotropic approximation, we shall consider only radially-symmetric potentials V(r) = 
V(|r|) for simplicity and their major occurence in physics. The trial frequencies decompose naturally into a radial 
frequency £Il and a transverse one SIt (see Ref. pj): 



! V — L ^2 h T I °MW 



(7.4) 



with r a = |r Q |. For practical reasons we rotate the coordinate system by x„ = C/x„ so that r a points along the first 
coordinate axis, 



(ra) 



r a , M =1 - 
0, 2 < /i < D, 



and r2 2 -matrix is diagonal: 



n 2 = 



tt 2 L 



Q| 



Qf, 



= un 2 u~ x . 



y ••• 

After this rotation, the anisotropic n-th order smearing formula in D dimensions reads 

i-h/3 



\ "^int [ r ] , 



( 27r )-^(«+i)/2(det al)- 1 / 2 (dot a 2 ,)-^- 1 )/ 2 



exp |~4 E 



z lk a Lkl Z U 



I 1 D n 

ex pi-oE E 



fj,=2 k,l = l 



The components of the longitudinal and transversal matrices a\ and a\ are 



(7.5) 



(7.6) 



(7.7) 



a Lki = a K T k,Ti), a% kl = a?r(T k ,Ti) 



(7.; 



where the frequency f2 in (5.13) must be substituted by the new variational parameters Ql, ^t, respectively. For the 



harmonic density in the rotated system we find 



Po L - T (r) 



-D-l 



27Ta loO V 2lTa T00 



exp 



2 a 



LOO 



2 a 



1 D 

2 M 
TOO „ =2 



(7.9) 



which is used to normalize (7.7) 
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The anisotropic smearing formula ( |7.7| ) will be applied to the Coulomb problem below. The anisotropy becomes 
significant only at low temperatures, where radial and transversal quantum fluctuations have quite different weights. 
The effect of anisotropy disappears completely in the classical limit. 



VIII. APPLICATIONS 

In this section we apply the theory to calculate the electron density of the hydrogen atom. For simplicity, we shall 
employ natural units with K = ks = M — 1. In order to develop some feeling how the approximations work, we first 
determine the particle density in a double- well potential to second order. 



A. The Double- Well 



In the case of the double-well potential 



(8.1) 



with coupling constant g, we obtain for the expectation of the interaction (6.3) to first order, also setting uS 2 = 1, 
< A„tN = i/3 5o + \giC$ ) H 1 ((x a - x m )/^2^j + \g 2 cf ] H 2 Ux a - x m )/^2^\ 



~8 93C P )H3 ~ x ™)/\f^A + J^94C { 4) H 4 ( (x a - x m )/\l2al Q 



•2) 



with 



3 111 

.90 = -«oo(^ 2 + 1) + 2-9 a oo + Sgaoo^m + y? x m + 2g~ ~ 2^™ 

.91 = -yfeaftoXm + ^g{2a 2 00 f /2 x m + g^j2^ x 3 m 
.92 = -aoo(^ 2 + 1) + 3 .9 a oo + SffOoo^m 
93 = g(2alo) 3/2 x m 



94 = 9%o- 



Inserting (3.2) in (|6.9|), we obtain the unnormalized double- well density 



Pi ' {Xa) 



exp[-/3Wi' Xm (x a ) 



(8.3) 



with the first-order effective classical potential 
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W ± ' (x a ) = -\n — +^{x a -x m ) t&nh— + - (A int [x}) XaXa . (8.4) 

After optimizing W^' Xm (x a ), the normalized first-order particle density pi(x a ) is found by dividing pi(x a ) by the 
first-order partition function 

Zi = -=^= / dx Q exp[-/3Wi(a; a )]. (8.5) 



Subjecting Wj ' Xm (x a ) to the extremality conditions (4. IS ), we obtain optimal values for fi (% a ) and x m (x a ). Usually 



there is a unique minimum, but sometimes this does not exist and a turning point or a vanishing higher derivative 
must be used for optimization. Fortunately, the first case is often realized. Fig. ||| shows the dependence of the 
first-order effective classical potential W 1 ' Xm (x a ) at 3 = 10 and g = 0.4 for three fixed values of position x a as a 
function of the variational parameters fl 2 (x a ) and x m (x a ) in a three-dimensional plot and its corresponding density 
plot. Thereby in both representations, the darker the region the smaller the value of W 1 ' Xm . We can distinguish 
between deep valleys (darkgray), in which the global minimum resides, and hills (lightgray). After having determined 



roughly the area around the expected minimum, one solves numerically the extremality conditions (4-13) with some 
nearby starting values, to find the exact locations of the minimum. The example in Fig. ^ gives an impression of the 
general features of the minimization process. First we note that for symmetry reasons, 

and 

n 2 {x a ) = Q 2 (-x a ). (8.7) 
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b) X a 




FIG. 4. Plots of the first-order approximation W 1 ' Xm (x a ) to effective classical potential as a function of the two variational 
parameters Q. 2 (x a ), Xm(x a ) at g = 0.4 and (3 = 10 for two different values of x a - 



Some first-order approximations to the effective classical potential W\ (x a ) are shown in Fig. |5| obtained by optimizing 
in fl 2 (x a ) and x m {x a ). The sharp maximum ocurring for weak-coupling is a consequence of the reflection property 



(3.6) enforcing a vanishing x m (x a = 0). In the strong-coupling regime, on the other hand, where x m {x a — 0) ~ 0, the 



sharp top is absent. This behaviour is illustrated in the right-hand parts of Figs. ^| and[?] at different temperatures. 
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FIG. 5. 



x a at p — 




n 2 (x a ) 



a) 




0.00 



b) 



-0.01 




VV9 



FIG. 7. a) Trial frequency Q 2 (x a ) at different temperatures and coupling strength g = 10. b) Minimum of trial oscillator 
Xm(x a ) at different temperatures and coupling g — 10. 

The influence of the center parameter x m dimimiishes for for increasing values of g and decreasing height l/4g of 
the central barrier. The same thing is true at high temperatures and large values of x a , where the precise knowledge 
of the optimal value of x m is irrelevant. In these limits, the particle density can be determined without optimizing in 



x m , setting simply x m = 0, where the expectation value Eq. (3.2) reduces to 



\ -^int 



\cf ] H 2 Ua/y^Q (ffl + 3. 92 ) + 92 cfHt Ua/fi^) + P (j9l + |<?2 + 9s) , (8.8) 



with the abbreviations 



9i = -«oo(^ 2 + !)) 92= gao Ql g 3 = —■ 



(8.9) 



Inserting (3.8) in Q6.9| ) we obtain the unnormalized double- well density 

1 



with the first-order effective classical potential 



/2^3 



expl-VWl'ixa)] 



(8.10) 



, . i sinh/?^ n 9 an 1 , . , lin 

(i.) = - In — ^- + -x 2 a tanh tL. + - ( ^ lnt [*] 



(8.H) 



The optimization at x m = gives reasonable results for moderate temperatures at couplings as low as g = 0.4, as 
shown in Fig. || by a comparison with the exact density obtained from numerical solutions of Schrodinger equation. 
An additional optimization in x m cannot be distinguished on the plot. An example where the second variational 
parameter x m does become important is shown in Fig. 0, where e we compare the first-order approximation with one 
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(Q) and two variational parameters (CI, x m ) with the exact density for different temperatures at the smaller coupling 
strength g = 0.1. In Fig. ^| we see that for x a > 0, the optimal x m -values lie close to the right hand minimum 
of the double-well potential, which we only want to consider here. The minimum is located at l/^/g « 3.16. We 
observe, that with two variational parameters the first-order approximation is nearly exact for all temperatures, in 
contrast to the results with only one variational parameter at low temperatures (see the curve for (3 = 20). Also 
for a small valuein the case g = 0.4, the optimization in fl 2 only gives reasonable results in the temperature region 
away from high- and low-temperature limits, as shown in Fig. |8| in comparison with the exact density obtained from 
numerical solution of Schrodinger equation. An optimization in two variational parameters gives no better result. A 
very instructive example, where the second variational parameter x m becomes important is shown in Fig. ^. There, 
we compare the first-order approximation for the double-well density with one (£1) and two variational parameters 
(f2, x m ) with the exact density for different temperatures at a coupling strength g = 0.1. This value of g no longer 
allows to neglect x m as in the case g = 0.4. We see from Fig. ^ that for x a > 0, the optimal x m -values lie close to the 
right hand minimum of the double-well potential, which we only want to consider here. The minimum is located at 
1/^/g ~ 3.16. We observe, that with two variational parameters the first-order approximation is nearly exact for all 
temperatures, in contrast to the results with only one variational parameter at low temperatures (see the curve for 
/3 = 20). 



0.2 
0.1 
0.0 

0.0 1.0 x a 2.0 3.0 

FIG. 8. First-order approximation using different for ft = 10 and g = 0.4 compared with the exact particle density in a 
double- well from numerical evaluation of Schrodinger equation. All values are in natural units. 
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FIG. 9. First-order particle densities of the double- well for g — 0.1 obtained by optimizing in two variational parameters 
Q 2 ,x m (dashed curves) and with only Q 2 (dash-dotted) vs. exact distributions (solid) for different temperatures. The parameter 
x m is very important for low temperatures. 



In second-order variational perturbation theory, the differences between the optimization procedures using one or 
two variational parameters become less significant. Thus, we restrict ourselves to the optimization in Q(x a )- 
The second-order density 

1 



exp[-/3W^(x a )] 



with the second-order approximation of the effective classical potential 
W?{x a ) = -\n 



i sinh/3n n 2 (m 1 n _ J_ , A? r n « 



(8.12) 



(8.13) 



requires evaluating the smearing formula (5.10) for n = 1 which is given in (S.8) and n — 2 to be calculated. Going 



immediately to the cumulant we have 
hp hp ( 

1 



<4 ? „tM ) XaW = dn dr 2 -(Q 2 + l) 2 [I 22 (n,r 2 ) - I 2 (n)I 2 (r 2 )} - -g(Q 2 + 1) [/^(n,^) - I 2 (n)h(r 2 )} 



+t^9 2 [hi(n,r 2 ) - h(n)h{T 2 )\ 
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with 



ImiVk) = (aoo - a 0kY 



exp 



j 2 + 2x a al k j 
On 2 (n 4 - a 4 ) 



3=0 



k = 1,2 



(8.14) 



(8.15) 



24 



and 



imn(n,r 2 ) = (-detA) 



dj™ dj 



n ex P 



2ag (detyl) 2 



(8.16) 



31=32=0 



det A = 



1 9 2 2 2 2 ('4,4,4n 
' Za 01 a 02 a 12 — a 00l a 01 "+" a 02 + a 12>- 



The generating function is 

F(ji,h) = OqoOi + jl) - 2al Q (al 1 j 1 + al 2 j 2 )x a + 2al {a\ 2 j 1 j 2 + (a^ + a^ 2 + af 2 )( a oiJi + Ooa^afa) 

/„2 • 1 2 • \/ 2 • 1 2 ■ 1 a 2 2 2 \ 



.17) 



P2(Xa) 




FIG. 10. Second-order particle density (dashed) compared with exact results from numerical solution of Schrodinger equation 
(solid) in a double-well at different inverse temperatures. The coupling strength is g = 0.4. 



All necessary derivatives and the imaginary time integrations in (S.14) have been calculated analytically. After 



optimizing the unnormalized second-order density ( 8.12 ) in f2 we obtain the results depicted in Fig. IfJ. Comparing 
the second-order results with the exact densities obtained from numerical solutions of the Schrodinger equation, we 
see that the deviations are strongest in the region of intermediate (3, as expected. Quantum mechanical limits are 
reproduced very well, classical limits exactly. 
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B. Distribution Function for the Electron in a Hydrogen Atom 



With the insights gained in the last section we are prepared to apply our method to the more physical problem of 
an electron in a hydrogen atom, with the attractive Coulomb interaction 



c 2 



V(r) = . (8.18) 



T 



Apart from the physical significance, the theoretical interest in this problem originates from the non-polynomial bature 
of the interaction. This makes the above-developed smearing formula is essential for finding variational perturbation 
expansions. Restricting our attention to the first-order approximation for the unnormalized density, we must calculate 
the harmonic expectation value of the action 

hp 

A„tM - J dn^t(r(ri)) (8.19) 
o 

with the interaction potential 

M„t(r) = - (^r + ^r T n 2 r\, (8.20) 



where the matrix has the form (7.4). We do not consider three more variational parameters r m here, because this 
will not be relevant in a strong-coupling case like the Coulomb interaction, as we know from the last section. After 
optimization in Q, we compare our results for the radial distribution function 



g(r) = p(r) (8.21) 

with the precise numerical results of Storer JlJ] . 

For the Coulomb potential, the optimization procedure can be simplified by setting the second optimization pa- 
rameter x m equal to zero from the outset. This is justified by observation made for the double-well potential, that 
the importance of knowing x m diminuishcs for decreasing height of the central barrier. Since the Coulomb potential 
has no central barrier, we may set x m = 0. 

1. Isotropic First-Order Approximation 



Applying the isotropic smearing formula (fTjf) for N = 1 to the harmonic term in (3.19) we easily find 



( r2(ri) )^ r =3 ^L_< + 4i r 2. (8 . 22) 
" a °00 °00 



2(> 



For the Coulomb potential we obtain the local average 

-^erf 



r(n) 



'oi 



r a a Q1 \ v/2aoo( a oo ~ «oi) 



.23) 



The time integration in (B.19) cannot be done in an analytical manner and must be performed numerically. Alterna- 



tively we can use the expansion method introduced in Subsection VI A for evaluating the smearing formula in three 
dimensions which yields 



/ a r i\" r at M-i e ra/2a °° H 2n+1 (r a / y/2a% ) r 2n ) f , T/ , [7TY \ - y 2 v , ^ 
- 4 -l r r, = yvr(J-,l] -^T^- 1^ 2 2n+i( 2n+1 )! °P J dy y V iat (^J 2a 2 00 y)e * H 2n+l {y). 



.24) 



o 



This can be rewritten in terms of Laguerre polynomials L%(r) as 



;Antw; 



• oo 
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FIG. 11. Radial distribution function for an electron-proton pair. The first-order results obtained with isotropic (dashed 
curves) and anisotropic (solid) variational perturbation theory are compared with Storer's numerical results (dotted) and 
an earlier approximation derived from the variational effective potential method to first order in Ref. Q (dash-dotted). 

Using the integral formula @ Eq. 2.19.14.15] 

OO 

dxx^e-^LKcx^cx) = (1 + 7)m(A ~ , Q + 3 F 2 {-m, a, a - A; 7 + 1, a - A - n; 1), (8.26) 

m\n\c a 
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where the (a) n are Pochhammer symbols, p F q (a\, . . . , a p ; b\, . . . ,b q ; x) denotes the confluent hypergeometric function, 



and r(x) is the Gamma function, we apply the smearing formula to the interaction potential (8.20) and find 



, i .<< e- x - (-l)"(2n-l)!! (2 „) , , , 

^nt[r]) ro , ro =-^rg 2 "(2n + l)! ^ H ^ + xK/ ^a^) 



n=0 

3 m~Z t 1 f „roi „ , , /„ n \ . 1 



- iV '2a6 m_ jc^^Cr./^ago) + -C^^/y^) j . (8.27) 
The first term comes from the Coulomb potential, the second from the harmonic potential. The resulting first-order 



isotropic form of the radial distribution function (8.21), which can be written as 



ff P(r„) = exp[-/?Wf (r )] (8.28) 

with the isotropic first-order approximation of the effective classical potential 

o , „ 3 sinh Btt O BQ, 1 , . . , , o , 

W?(r a ) = —In + -rl tanh P — + - (A int [v] )£ ^ , (8.29) 

is shown in Fig. [H] for various temperatures. The results compare well with Storer's curves |]l4j| . Near the origin, 
our results are better than those obtained with an earlier approximation derived from lowest-order effective classical 
potential W\{xq) Q. 



2. Anisotropic First-Order Approximation 



The above results can be improved by taking care of the anisotropy of the problem. For the harmonic part of the 



action (3.19), 



A int [r] =An[r}+A c [r] 



(8.30) 



the smearing formula (7.7) yields the expectation value 



\ {^i«ioo [cf + ^S^Cra/v^)) + 2^4 00 (c( 0) - c«)} , (8 



3D 



where the C^ L , T ^. are the polynomials ( |6.5|) in which Q is replaced by the longitudinal or transverse frequency. For 



the Coulomb part of action, the smearing formula (7.7) leads to a double integral 



h/3 



[Mr] 



7ra loo( 1 - a i) 



a !oo( 1 - a i) 



exp 



.32) 



with the abbreviations 
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2 _ LQ1 2 _ "TQl f o oo\ 

a L — — — , a T — — — . (6.66) 

a L00 a T00 

The integrals must be done numerically and the first-order approximation of the radial distribution function can be 
expressed by 

g^ T (r a )=eKp[-pW^- T (r a )] (8.34) 

with 

This is optimized in r2i(r a ), J7T(r a ) with the results shown in Fig. [Oj The anisotropic approach improves the isotropic 
result for temperatures below 10 4 K. 



IX. SUMMARY 

We have presented variational perturbation theory for density matrices. A generalized smearing formula which 
accounts for the effects of quantum fluctuations was essential for the treatment of nonpolynomial interactions. We 
applied the theory to calculate the particle density in a double-well potential, and the electron density in a Coulomb 
potential, the latter as an example for nonpolynomial application. In both cases, the approximations were satisfactory. 
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